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Abstract 

We address the problem of the so-called "granular gases" , i.e. gases of mas- 
sive particles in rapid movement undergoing inelastic collisions. We introduce 
a class of models of driven granular gases for which the stationary state is the 
result of the balance between the dissipation and the random forces which 
inject energies. These models exhibit a genuine thermodynamic limit, i.e. at 
fixed density the mean values of kinetic energy and dissipated energy per par- 
ticle are independent of the number N of particles, for large values of N. One 
has two regimes: when the typical relaxation time r of the driving Brownian 
process is small compared with the mean collision time r c the spatial density 
is nearly homogeneous and the velocity probability distribution is gaussian. 
In the opposite limit r 3> r c one has strong spatial clustering, with a frac- 
tal distribution of particles, and the velocity probability distribution strongly 
deviates from the gaussian one. Simulations performed in one and two di- 
mensions under the Stosszahlansatz Boltzmann approximation confirm the 
scenario. Furthermore we analyze the instabilities bringing to the spatial and 
the velocity clusterization. Firstly, in the framework of a mean-field model, 
we explain how the existence of the inelasticity can bring to a spatial clusteri- 
zation; on the other side we discuss, in the framework of a Langevin dynamics 
treating the collisions in a mean-field way, how a non-gaussian distribution of 
velocity can arise. The comparison between the numerical and the analytical 
results exhibits an excellent agreement. 

PACS: 81.05.Rm, 05.20.Dd, 05.40. +j 

I. INTRODUCTION: HYDRODYNAMICS AND GASES 

Granular systems (sand, powders, seeds, cements, etc.) have been extensively studied, 
in the last two decades, by means of analytical investigations, experiments and computer 
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simulations. The rich and intriguing phenomenology is well known to engineers and the need 
of a better comprehension of granular behaviors is widely recognized in applied sciences as 
well as in theoretical physics. A quite exhaustive review may be found in M. Problems 
in granular systems are roughly divided in quasi-static (sand piles, distribution of static 
forces, compaction, fractures propagation, etc.) and dynamical ones (all kind of flows, 
convection and segregation, pattern formation, ffuidized beds, etc.). In the latter class, 
large collections of inelastic particles are involved in fluid-like rapid dynamics, therefore the 
hydrodynamics approach seems to be the natural one. The main granular hydrodynamics 
theories are reviewed in [jj: they are all based on non-equilibrium conservation laws || for 
mass, momentum and energy: 

t + pV-u = (1) 



p _ = _ V .r + « (2) 

l,| = _ V . q + f:Vu-r (3) 

where p is the local density, u is the local velocity vector, d/dt = d/dt + u ■ V is the La- 
grangian derivative, f is the total stress tensor (f = f s + f c , as transport of momentum has 
two contributions: a "streaming" term and a "collisional" one), g is the body-force vector 
(gravity or else), T is the so-called granular temperature T =< (u— < u >) 2 >, q is the 
flux vector of "granular heat" , f : Vu is the tensorial product for the granular-temperature 
generation by shear work and F is the sink term due to dissipation into thermodynamic heat 
(i.e., energy lost in collisions). The existing approaches differ in the constitutive relations 
that make r, q and F depend on the other properties Vu, p, T and on all the parameters 
of the system. Apart some heuristic relations (see, for example, Q or ||), there are many 
calculations based on the kinetic theory of nonuniform dense gases ||, making some assump- 
tion on the form of the single particle distribution function /(x, v, t), that is the solution 
of the Boltzmann-Enskog kinetic equation. Before j7|, and the almost simultaneous M, all 
the theories had assumed a Maxwellian velocity distribution, but the nonzero off diagonal 
components in the streaming stress tensor (strongly apparent at low solid fractions) indi- 
cate the need of a correction to Maxwellian distribution. In and || accurate predictions 
of streaming stresses, in the case of slightly inelastic and slightly rough smooth spheres, 
are obtained. In and || and in successive studies (see 0), a certain degree of energy- 
equipartition breaking is assumed, considering two different temperatures for translational 
and rotational degrees of freedom, respectively. 

The increasing power of computers awakened an interest in "granular gases" simulations, 
that is the investigation of kinetics of granular systems far from close-packing, e.g. a granular 
phase diagram has been proposed in |J to clarify different regimes and distinguish between 
a gas-like phase and two different high-density phases. Results from these simulations have 
been compared to previous granular hydrodynamics predictions, showing a disagreement in 
various aspects. 
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Many simulations have been performed in the "cooling" situation ( flOf , UTTf , [12], [13 

0, El 



the particles evolve with no external forcing, dissipating in collisions all initial 
energy. Hydrodynamic predicts a time-scaling of granular temperature T ~ t~ 2 under the 
assumption of Maxwellian velocities at all times, but strong departures from this law are 
observed when (at fixed volume) the number of particles grows or when the restitution 
coefficient r decreases ||11|| . When N(l — r) >> 1 |L2) (where N is the number of particles) 
it is found that the clustering instability (this can be derived [TI| from Jenkins & Richman 
hydrodynamics ||) may degenerate in the so-called inelastic collapse, as particles may be 
trapped in a sequence of infinite collisions in a finite time (i.e. a divergence of collision 
rate). 



13 1 simulations. 



Inelastic collapse is found in one-dimensional and two-dimensional 
Furthermore, equipartition between rotational and translational energy is found to be broken 
in cooling kinetics (see |0| and references therein). 

Strong equipartition breaking is found in another class of models, that of driven granular 
gases, where the dissipation of energy due to collisions is balanced by an external source 
(in realistic situations one has to vibrate or shake a granular system, to keep it alive). 
We could divide these models in two sub-classes: elitary and democratic models, referring 
to the quantity of particles receiving energy from the external source. In elitary (one di- 
mensional |TB[, [1711 and two-dimensional [|T7|) models a wall of the container is the unique 
energy source, therefore there are few particles (just one, in Id) that transport energy from 
the source to the rest of the system. In democratic models (one dimensional [pl|, [l9j and 
two-dimensional [[2(J) all particles receive energy, by mean of a Brownian-like random ve- 
locity kick at every time step. The model we propose in this paper belongs to this last 
sub- class. 

Recently the Boltzmann-Enskog equation for granular kinetic (cooling or driven) has been 
analyzed, showing that the velocity distribution is expected not to be Maxwellian. Esipov 
and Poschel [^TJ have found, for cooling inelastic hard spheres, exponential tails while van 
Noije and Ernst [G2J have obtained the same tails in the cooling regime and ~ exp(— Av 3 ^ 2 ) 
tails in the (democratic) driven regime. 

A non-Maxwellian behavior has been, very recently, observed experimentally in a verti- 
cally driven granular bed [ 2"3"| . The measured velocity distribution in such experiment is in 



very good agreement with the results of our simulations. 

To conclude this brief introduction, one has to remember that, in modeling granu- 
lar gases, the oversimplifying criterion is, sometimes, misleading. Brilliantov et al. have 
shown \ Z IA\ that the universally accepted picture of fixed restitution coefficient is far from 
being obvious and that the behavior of granular gases (self diffusion, as an example) may 
change drastically if this coefficient is taken dependent on the impact-velocity. McNamara 
and Luding have recently stressed the importance of keeping into account rotational de- 
grees of freedom as well translational ones, and that of the particle-wall dissipation. Further- 
more, the relevance of choosing the boundary energy source and the possibility of considering 
soft particles instead of hard grains have been investigated by Geisshirt et al [^6 . 



In section [TT] detailed results of simulations of a one-dimensional model of driven granular 
gas (already presented in p7| ) are reported. In section [III] a one-dimensional and two- 



dimensional version (in Boltzmann approximation) of the same model is discussed, reporting 
analogous results. Then, in sections [TV] and [V], some theoretical interpretations are proposed 
in order to understand non-maxwellian behavior and try to relate it to the clusterization 
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phenomenon. Section [VI] is devoted to conclusions and open problems. 



II. THE INELASTIC HARD-ROD ONE-DIMENSIONAL MODEL 

Such class of models originates from the seminal paper of Du, Li and Kadanoff |nj who 



considered N identical hard rods confined between a thermal and a reflecting wall. In this 
model one has a statistically steady state as the result of the balancing between the dissipa- 
tion of the kinetic energy due to the collisions between the rods and the energy reinjection 
due to the thermal wall; the latter supplies energy only to the last particle, which in turn 
transfers energy and momentum to the rest of the system, producing a somehow trivial 
cluster near the opposite wall. Such a state represents a breakdown of the equipartition of 
the energy in a stationary non-equilibrium system; however, its existence comes about as an 
artifact, since it is due to the peculiarity of the boundary conditions. In fact, in the model 
introduced in ref. |TB| the mean kinetic energy per particle, 



TV 



8=1 

and the mean dissipated energy per particle per unit time, 
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(where (Ai£)j is the energy loss during the j t h collision occurred in the time interval [t — 
At/2,t + At/2], and < > is the time average), decay exponentially with the number N of 
particles, as shown in figure ([I]). Finally, within the Kadanoff model only a small region of 
parameter space can be explored: since r, the restitution coefficient, defined below, must 
satisfy the inequality N(l — r) < 1 in order to avoid inelastic collapse flO |, [ |12[ . 



Williams et al. |18[ proposed an alternative heating mechanism. The idea is to supply 
kinetic energy to every particle by means of a random acceleration at every time step. Since 
the dissipation due to inelastic collisions is not effective in balancing the increase of energy 
coming from the random kicks (the latter is independent from the velocities, while the 
former is proportional to them), the authors subtract the average velocity of the center of 
mass of the system from the velocity of each particle at every time step in order to avoid 
energy divergence and total non-conservation of the total momentum. Though this method 
is numerically efficient, it does not appear realistic from a physical point of view. 

Hereafter, we propose a model |27j consisting of N identical particles of mass 1 on a ring 



of length L. Between collisions, each particle obeys to the following Langevin equation ||28|| : 



dVi Vi , 2Tf s(+\ fa\ 

-dF = -- + y— m (6) 



% = (7) 
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where 1 < i < N, r and Tp are the relaxation time due to viscous effects and the thermal 
bath temperature respectively, fi(t) is a standard white noise with zero average and < 
fi(t)f j {1f)>=S ij 6(t-tf). 

In addition to these equations, the particles mutually collide according to the following 
rules: 

a) only binary collisions are considered, 

b) each collision is instantaneous, 

c) the post-collisional velocities are related to the pre-collisional ones by the equations: 



, 1 — r 1 + r 

,1+7-1 — 7- /nN 

Vj = -^r v i + ~^r v o (8) 



where r is the restitution coefficient. In this way, the momentum is conserved in the 
collisions, while the kinetic energy of center of mass is rescaled by r 2 , i.e. 

(vl-vtf =r 2 ( Vi - Vj ) 2 (9) 

The elastic case is r = 1, while for r = the colliding particles have no relative motion after 
the collision: they move together with the velocity of the center of mass. It has to be noted 
that, in one dimension, the size of the particles is not a significant parameter, because of the 
"hard" nature of collisions: the particles never deform (this effect is kept into account in the 
restitution coefficient picture) and therefore only the length of spacings between particles is 
important, that is L (see fHf] ) 

In absence of collisions, each particle would perform a Brownian motion reaching, for 
t » r, a stationary state with a gaussian velocity distribution: 



or 



p ^=vm exp {-*k ] (10) 



and a diffusive behavior: 



<(x 4 (t)-x 4 (0)) 2 >=2 J Dt (11) 

where D = Tpr/2 is the diffusion coefficient. The interpretation of the model is straight- 
forward: the environment supplies kinetic energy to the system as a thermal bath at tem- 
perature Tp. The viscous term (with characteristic time r) is naturally introduced to take 
into account different friction effects, such as contact friction with boundaries, particle-fluid 
interaction, tangential inter-particle friction and energy transfer among different degrees of 
freedom. Experimental fluidized beds, see for example ||29||, are an example of phenomena 



in which a viscous damping and a noisy term are naturally present. As noted before, in [23 
is presented an experiment showing strong analogies with our model. 

When collisions are considered, another characteristic time emerges, that is the average 
collision time r c between two successive encounters. An estimate of function of 

average density and typical velocity, is 
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Tc ; (12) 

2NV< v 2 > 

It is natural to assume that < v 2 > reaches a stationary value with statistical fluctuations 
(of order ~ 1/iV), as it is observed in simulations. In the following we shall refer to the 
quantity 



N t 

r » =< » 2 >= L < 13 > 



T 



as to the granular temperature of the system. Note that the system is not at equilibrium, 
therefore T g is not a temperature in a proper thermodynamic sense: it may be different, if 
one measures it at different scales or in different subsystems (as it will be shown later). 

In all the simulations performed, we use L/N = 1 and T F — 1 and the measured T g < Tp 
is almost always found greater than 0.1 (and never drops below 0.03). ^From eq. (|I2| ) we 
can estimate 0.5 < t c < 5. 

The presence of two time scales (r and r c ) leads to two different regimes. As r c varies in 
a small range (less than one order of magnitude), we could tune the parameter r to observe 
these two phases: 

A) When r << r c the effect of collisions is rapidly overwhelmed by the Brownian motion 
(i.e. collisions are rare events) and the system behaves as a collection of weakly interacting 
random walkers or, equivalently, as an ideal gas in equilibrium at a temperature T g not far 
from the temperature Tp of thermal bath (one can be convinced of this also looking to the 
Boltzmann equation of the system, reported in the next section). 

B) In the opposite limit r c << r the collisions dominate the dynamics and strongly 
compete against the driving mechanism. In this regime a statistically stationary state is still 
observed, in the sense that macroscopic averages are well defined, and interesting phenomena 
emerge: 

a) a strongly inhomogeneous spatial distribution {clusterization) 

b) a deviation of velocity distribution from gaussian behavior. These phenomena are 
more and more pronounced with decreasing values of the restitution coefficient r. 

The simulations have been performed using a fixed step At integration of eqs. (||) and 
(0) where At << r c and an event driven check of collisions during every time step. For 
low values of r we observed an exponential decreasing of collision time, much shorter than 
the integration time At. We discarded these simulations, interpreting them as examples of 
inelastic collapse |12j| . The critical value of r, for the appearance of collapse, increases with 



r: in the limit r — > oo the thermal bath disappears and the system becomes a granular 
cooling model with critical value, for inelastic collapse, estimated by iV(l — r) ~ 1 as noted 
before. 

In figure (Q) we report T g and W vs. the restitution coefficient r for different r. 
A simple relation between T g and W may be obtained. The variation of the kinetic 
energy due to Langevin dynamics is 
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N N 

(SE(t)) Lang = — ^(t) + 5 Vl (t)Y - — = (14) 

i=l i=l 

1 - 1 - 

i=l i=l 

where 5vi is the velocity variation during a time interval dt in equation @, from which we 
obtain the relations: 



|u x ' v i {t)s Vi {t) \_ ({vm?) (16) 



dt-»0 \ dt 

where the < ... > average is taken over different realizations of stochastic process fi(t). 
Recalling the definition of T g , using equations (||), (|5[), and inserting eqs. (|T|) and ( |TT)| ) into 



eq. (|T4]), and assuming the ergodicity, one obtains: 



r 

The numerical check of such relation is shown in figure (§) 

Though the system is statistically stationary, the instantaneous density of particles is 
rapidly evolving. To get an idea of different density profiles in the two regimes (homogeneous 
and clusterized) , look at fig.®. 

The density distribution becomes fractal in the clusterized regime, as can be verified 
measuring the correlation dimension d 2 @] that we computed from the correlation function 



C(R) = N2 ^_ to) j t dtJ2v(R-\x*(t)-Xj(t)\)~R d2 (18) 



with to the time after which one can assume the system is in a typical situation and T is 
the duration of the simulation. C(R) is shown for homogeneous and clusterized situations in 
figure (|5[). In figure (Q) is presented a summary of d 2 measurements as function of parameters 
t and r. 

The clusterization may quantitatively characterized by means of an entropy defined as 

M 
3=1 

where the ring of length L is divided in M equal boxes (i.e., segments) and rnj is the 
number of particles in the j t h box. The entropy hu attains its maximum value hu = logM 
when rrij = N/M for every box j. hu decreases as the density distributions becomes more 
and more clusterized. For a non-clusterized (but fluctuating) density, we have a Poisson 
distribution for rrij, that is (with A = N/M) 

/K) = — exp(-A) (20) 
rrij. 
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from which it can be (numerically) calculated the effective entropy for homogeneous regime 
h* M . In figure (0) are presented many measurements of H M jH* M where H M = exp(< h M >), 
H* M = exp(< h* M >) and < > is the time average. The quantity H M / H* M basically gives an 
indication of the fraction of non-empty boxes in a typical snapshot. 

In figure @ is shown the distribution of velocities, obtained sampling the velocities of 
all particles for very long times, in the two different regimes (a quasi-equilibrium case with 
t = 0.01, r = 0.99 and an out of equilibrium case with r = 100, r = 0.7). In the quasi- 
equilibrium regime the distribution is very well fitted by a gaussian. As a general result, 
when t >> t c the velocity distribution ceases to be gaussian and the high velocity tails decay 
more slowly to zero. The deviation becomes more pronounced as the restitution coefficient 
r decreases. In the figure a theoretical fit is also plotted for the non-gaussian distribution. 
This fit will be discussed below, in section [TV[ 

In figure @ the density distributions in the homogeneous and clusterized regimes are 
shown (respectively r = 0.01, r = 0.99 and r = 100, r = 0.7, with iV = 300), i.e. fuijn) 
where m is the number of particles in a box when the ring is divided in M boxes (this 
distribution, as all the others, is obtained sampling data for very long times). The homo- 
geneous regime is very well fitted by a Poisson distribution, as noted in equation (|2Cf) . The 
clusterized regime presents an exponential long tail and a power law for the low density 
boxes: the function -Le~ cm with c = 0.14 fits very well the data and is consistent with the 

m J 

theoretical interpretation given in section |V[ 

The figure (|10D represents the box granular temperature T/vf(m) as a function of the 
number of particles m in a box: 

j m 

T M (m) = - X>;- < v > m ) 2 (21) 

m 3=1 

where < v > m is the average velocity (typically close to zero when m » 1) in the box and 
M is the number of boxes on the ring of length L. The figure shows this function (averaged 
over very long times) for both the regimes: in the gaussian case (with r = 0.01 and r = 0.99) 
we observe that T M {m) is a constant, while in the non-gaussian case (r = 100, r = 0.7) it is 
a power law, i.e. Tm(tti) ~ m _/3 with (3 = 0.5. The exponent (3 depends on the values of r 
and r. 



III. THE BOLTZMANN EQUATION APPROXIMATION 



A natural question now arises: can we expect that the above results are general and 
independent from the dimensionality of the system? Or are these an artifact of the one 
dimensional dynamics? An answer may come from the Boltzmann equation for the one 
particle distribution P(x,v,t) plf : 

1 d(vP) T F d 2 P _ dP 
= ~dt 



BP d{vP) 
dt dx 



dv 



dv 2 



coll 



dP 
~dt 



4A 



coll 



dv'\v' - v\P(x, v\ t)P(x, (2v - (1 - r)v')/{l +r),t)- 
— A / dv'P(x, v ', t)P{x, v , t) \v' — v | 



(22) 



(23) 
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where A ~ l/r c is the mean collision rate per particle. In the limit of elastic collisions 
(r = 1) the collision integral (|23| ) disappears and the stationary solution of eq. ( ^2|) is 
P(x,v,t) oc exp(— v 2 /2Tp). This is related to the fact that in the elastic limit a collision 
between two particles is nothing but a change of the labels of the two particles and therefore 
the collisions are not relevant at all. 

The main approximation in eq. (|22|) is the Boltzmann Stosszahlansatz, according to which 
the correlation between two close particles is neglected and one writes: 

P 2 (x, x', v, v', t) = P(x, v, t)P(x', v', t) (24) 

As equation (p2|), as far as we know, cannot be solved analytically, we consider a stochas- 
tic process based on the Bird algorithm ||32|| , whose statistical features are identical to those 



of the Boltzmann equation. N particles move on a torus (in d dimensions) of linear size L 
(i.e., the area of the d-torus is L d ). The time is discretized in intervals of duration At. A 
collision time r c is fixed a priori: this means that two particles collide, during At, with a 
probability p = At/r c (At << r c , as usual). 

At each discrete time t^ = kAt, positions and velocities are upgraded according to eqs. 
(^) and (0). Then for each particle i a random number y is extracted out of a uniform 
distribution in the interval [0,1]: if y > p no collision occurs, otherwise the particle i 
collides with a particle j such that \xi(tk) — Xj{tk)\ < I (I << L), chosen with probability 
proportional to |fi(tfc) — Vj(tk)\. The collision rule, in one dimension, is the same as before 
and it is extended to the two and three dimensional cases in a natural way: after a collision 
Vi'— vj'= re(vi— vj) where e is a unit vector with random orientation. 

It can be demonstrated JHJ that for this process, in the limit N — > oo, p — > 0, I — > 0, 
At — > 0, the evolution of the probability distribution P(x, v, t) of the Bird process is governed 
by the Boltzmann equation ([22]) . 

In this model there are two parameters, r c and I. The first one was already an observable 
of the previous model (where it was almost constant). The second one, instead, represents 
the collision range and replaces the radius of the particles, so it is related, in some way, 
to the total volume fraction v = kN/L d with k the volume of one particle, which is not 
explicitly considered in this model and in the previous one. Furthermore, we cannot expect 
to find a fractal scaling at a range lower than I, and, moreover, no inelastic collapse can 
be observed, as short range correlations are neglected. It should be noted, finally, that the 
imposed collision time r c is larger (but of the same order of magnitude) than the one really 
measured in the simulations, r*, because a particle collides with probability p = At/r c only 
if there is another particle at distance lower than /. One can expect that r* — > r c as the 
clusterization becomes more and more pronounced. 

In all the simulations of this model we kept Tp = 1 and N/L d = 1, we rarely changed r c 
and / and we explored the space of parameters r and r as in the previous section. All the 
results previously obtained were recovered in this approximation with 1 < d < 3, showing 
that they are general properties of a granular system subject to such a heating mechanism 
as that of eqs. (0). A rapid overview of significant measurements, in one and two 

dimensions, follows. 

We stress the fact that also in this model the system reaches a statistically stationary 
state after a transient, and also in this model there are two different regimes: the quasi- 
equilibrium regime (r << r c ) and the out-of-equilibrium one (r >> r c ). 
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The clusterization phenomenon is shown in figures (|TT|), flT2|) and ( |13D where the density 
snapshots and the correlation functions C(R) (defined in eq. (|l8l)) are presented. It is 
observed the predicted reduction of the fractal scale (more evident in the two-dimensional 
model) due to Boltzmann approximation. 



The existence of a good thermodynamic limit is shown in the figures (|nj) and ( p~5j ) where 
the iV-dependence of some observables is plotted: we show the granular temperature T g , the 
fractal dimension 0*2 and the collision rate \ = 1/r*, all in the out-of-equilibrium regime. 

The same quantities are plotted in figures ([11]) and ([[7]) against the restitution coefficient 
r in the same regime (these plots are analogous to those of figures (|2|) and (|6|)). Note that 
t* approaches r c when r decreases, as it is expected. The diffusion coefficients, also plotted 
in those figures, will be discussed in the conclusion. 

The distributions of velocities are presented in figure ( |T8D and (|T9|). The non-gaussian 
behavior in the second regime (r >> r c ) is still clearly observed. In figure (|20"D a distribution 
of velocities restricted in the following way is presented: we sampled the velocities of one 
particle only when there were other m (m = 1 and m = 5) particles in a box of radius R 
(I << R << L) centered on that particle. This is a sort of velocity distribution at fixed 
density. The plot shows a less pronounced deviation from the gaussian, in agreement with 
the statement (discussed in section |TV| ) that there is a local equilibrium with a temperature 
that depends upon the local density, in order to have a stationary distribution of clusters. 

The analogue of the figure (H) and (|10D are the figures (|2l|) , ( f22|) , ( ^3|) and (|4]). Again 
the density distribution fu ijn) is a Poisson function when the system is not clusterized and a 
function of the form exp(—am)/m^ in presence of clustering. The box granular temperature 
Tjvf(m) still presents a constant behavior (as expected) at equilibrium and a power law m~ 7 
in the non-gaussian regime. It does not seem possible to easily find a relation between a, (3, 
7 and the other observables (as d 2 or T g ). 

In summary, the exact model is perfectly reproduced in this Boltzmann approximation, 
at large N, not only confirming the existence of a thermodynamic limit, but also showing 
that this system may be studied analytically in order to better understand this kind of 
driven granular kinetics. In the following sections some theoretical interpretations will be 
discussed. 



IV. A MODEL TO EXPLAIN CLUSTERIZATION 

In this section we address the problem of the microscopic origin of the clusterization. In 
order to do that, we study a class of models in which the system is composed by M boxes 
and N particles in a mean-field model, assuming that the boxes have infinite connectivity. 
One starts with a certain configuration and let the system evolve with an exchange dynamics 
in which, at each time step, one particle moves from one box to another, both boxes being 
chosen randomly. The probability for each single exchange is model-dependent and it will 
be our tuning-parameter to scan the different phenomenologies. Our goal is to understand 
in a quantitative way how the microscopic dynamics affects the clustering properties of 
the system. In particular we shall try to recover the results, obtained in the framework 
of the models previously introduced, for the density distributions in the clusterized and 
homogeneous cases (see Figs. fl2ip . (|22|)). 
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The models are defined in terms of master equations for the probability P m of having a 
box with m particles, assigning transition rates for landing in a box with m particles Wi n (m) 
and for leaving a box with m particles W out (m). It must be 

W m (N) = W out (0) = (25) 

and the normalization conditions must be satisfied: 

OO OO jy oo 

Y,Pm = l Y. mPm= M = X ^TW in (m)MP m = l (26) 

m=0 m=0 m=0 

The general question is: what is the asymptotic stationary distribution for the average 
number of boxes with m particles, P(m)7 

The simplest case we can consider is the one in which each single movement is indepen- 
dent of the state of the departing and of the landing box. In this case there is no bias in the 
movements and W in (m) and W out (m) do not depend upon m: 

W in (m) = W out {m) = (27) 

and the general master equation reads 

M 2 = .P m _i(P m _i — — ) + 2_P m+ iP m _i + P m+ i(P m+ i — — ) + (28) 

~hfm+l(l — Pm ~ Pm+1 ~ Pm-l) + (1 — Pm ~ Pm-l ~ Pm+1 ~~ P<o)Pm—l + 
P'm—lPrn ^Pm\Pm ) PmPm+1 P-rni^ Pm—1 Pm Pm+l) 

-(1 - P rn - P rn+l - P m _! - P )P m for < m < N 

HP 1 
M ._l = fl(1 _ ft __ ) _ (1 _ fi _ J y fl| 

M2d -w=™->- w-h) (29) 

In the limit of M >> 1 one can neglect the ^ terms in the right hand side of eq. (p§[) 
and easily get the stationary solution = 0) 

P m = Ae- cm (30) 

with A = 1 — e~ c corresponding to the normalization condition P m = 1 and where c is 
a constant depending on iV and M: c = ln(l + v) with A = -jr. 

This result has to be compared with the probability fuijn) m the non-clusterized case 
of the previous sections. In order to do this it is necessary to recall that this result has been 
obtained with a small value of the number of boxes M. This means that one is very far from 
the limit M » 1 and this situation corresponds to a sort of coarse graining in the system 
in which each box (big box) is actually composed by a certain number of small boxes (whose 
number is such that M » 1). The problem can thus be formulated in the following way: 
given a system of N particles distributed in M smaU boxes with the distribution P m given by 
eq. (|30"D , what is the distribution P^ for the particles in a system of M big boxes each one 
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composed by R (R = M sma u/Mug) small boxes? The resulting distribution is easily written 

as 



* R 



P rn = J2U P ^ = A R e- cm F(m,R), (31) 
i=i 

where ^* indicates the sum on the {mi,...,m#} such that J2f=i m i = m i F(m,R) is the 
number of ways of distributing m particles in R boxes and it is given by [j34jj : 

F(m,fl)=( m + f- 1 ) (32) 

With the help of (^) and using the Stirling formula, the expression ([JT]) becomes (for 
R » N » 1) 

P * _ A R r -c m (m + R-l)\ „ fi (m + fl-l)> + fl-r ^ 

^~ Ae m\(R-l)\~ Ae m\{R-l) R ^ ~ 1 ' 

~ A R e - cm — = ^g-C^d+^r)-^-^) w 



m! m! 

-A* l A J 

m! 



It has been used the definition of R, the fact that c = ln(l+M sma ^/A^) and that M smaU /N » 
1. In the last passage A* = N/Mug has been introduced and has become e~ A *, as can be 
verified when A -1 = M sma u/N » 1. It has been shown, therefore, that the coarse grained 
version of the solution of fl28D is exactly the Poisson distribution found in the simulations, 
in the non-clusterized regime (see figs. (£5|), ( pTj ) and (|22|)). 



Let us consider now one case where the transition rates for the particle jumps depend 
on the contents of the departing and landing boxes. This corresponds to impose some sort 
of bias to the system that could well reproduce the situation one has in the clusterized cases 
due to the inelasticity. We consider in particular the following case, defined by the transition 
rates: 

wwo) = JL 

W in (™) = (l-P )- forO<m<iV (34) 



Tfl 

W out (m) = — for < m < N 



N 

N 

These transition rates, that satisfy the relations (p6|) , have the following interpretation. 
The probability to land on a box containing already m particles is proportional to the 
number of particles because this mimics the inelastic collision with a cluster of m particles. 
On the other hand the departure from a box containing already m particles has a probability 
proportional to m because the probability to select one particle in that particular box is 
proportional to m. 

Neglecting as usual the terms of the order of and after simplifications, the stationary 
master equations write: 
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P,„ + i(m + 1) + (1 - P )(m - l)P„_i - P mP m = (35) 

|p n _ i( i-p )^-p w( i--L)± = o 



The solution in this case is given by: 



P m = A-e~ am P = 1 - e" a (36) 



m 

with A = X(e a - 1) and X = §. A and a are related by an implicit equation obtained 
imposing the condition P m = 1, that in the limit iV — > oo becomes 

l-e- a - A*\n{l-e- a ) = 1 (37) 

In the clusterized case we expect the solution to be self-similar, in the sense that P m has the 
same behavior of P^, and the coarse graining previously performed should not change the 
solution (PB[), apart a rescaling of A and a. 

It must be noted that, as A must be finite, when N — > oo (and M is fixed) a has to 
go to zero, while a diverges when N/M goes to zero. It is natural to think to a as to the 
inverse of the characteristic "mass" of a cluster, that is the typical number of particles in it. 
In this sense the term exp(— am) acts as a finite-size cut-off for the self-similar distribution 
P m ~ 1/m. 

The solution (|36"D is in excellent agreement with the numerical results obtained in the 
previous sections. In particular in the case N = 300 M = 100 of the one- dimensional model 
of Sect. II one recovers the density distribution with the correct value of a ~ 0.14 (see fig. 

(I))- 

To get the other observed behaviors of density distribution P m ~ e am jmP (see figs. 
(PI) and flUD), it is enough to change the transition rates appearing in eqs. (j34]) into the 



following: 

W in (0) = ± (38) 

W in (m) = n(l - P Q )m p for < m < N (39) 

Wout(m) = nwP for < m < N (40) 

where \x is a normalizing constant: 

^(M^py) . (4i) 



Now, we can go a step further relating the clustering properties of the system to the ve- 
locity distribution. In order to do that we consider the following quantities: the distribution 
of boxes, fu{ m \ containing a given number m of particles and the velocity variance Tj^(m), 
in a box occupied by m particles. We consider first the non-clusterized case (r << r c and 
r ~ 1). Within this regime we find from the simulations that: 
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T^ as {m) ~ canst. (42) 

\rrip—m 

fM S (m) = ——. (43) 

By assuming in each box a gaussian velocity distribution with a constant variance T^ as (m) 
it turns out that the global velocity distribution p elas (y} is gaussian. Let us recall that 
the Poisson distribution is the one associated with a process of putting independently AiV 
particles into iV boxes. 

Let us turn to the non-elastic case. If r = 100 and r = 0.7, considering the occupied 
boxes (m > 0), we obtain from the simulations the following relations 

T™ l (m)~m-P (44) 

„— am 

flT\m) = ——, (45) 

with (3 ~ 0.5 and a ~ 0.14. Let us compute from these scalings the global velocity dis- 
tribution. Taking into account that the spatial probability distribution of the particles is 
/m(^) and assuming that their local velocity distribution is gaussian, but with a variance 
Tjw(m) ~ m _/3 which depends on the occupancy, we obtain, for the global velocity distribu- 
tion Pi ne i(v), and in the continuum limit: 

p mel (v)^Y. ei ~ vI ^ )e ~ am - ( 46 ) 

m=l 

We stress how the the distributions measured in the simulations are in very good agreement 
(see the dashed line in fig. (]8p) with the numerical computation of eq. (fiSD, which, in 
summary, has been obtained under only the following hypothesis: 

(i) non-Poissonian distribution for the box occupancy fui 111 ) oc e~ am /m; 

(ii) gaussian distribution of velocities in each box with a density-dependent variance 
?Af(m) cx m _/3 . 

The hypothesis about the scaling relation between the velocity variance (i.e. Tw(m)) and 
the local density, apart from being justified numerically, can be understood in the following 
way. The stationarity and the scale-invariance of the cluster distribution, implies a certain 
distribution of lifetimes for the clusters. In particular each cluster has a lifetime which is 
inversely proportional to its size. The scale-invariant cluster-size distribution thus implies 
a scale-invariant distribution for the lifetimes. The cluster lifetime is strictly related to the 
variance of the velocity distribution inside the cluster itself. In order to ensure the stability of 
a cluster in a stationary state we have to require that the velocities of the particles belonging 
to it are not too different, or equivalently that the variance of the distribution is smaller the 
higher the density. So, given a scale-invariant distribution of clusters one would expect a 
scale- invariant distribution of variances, that is Tm(wi) ~ m" 13 . 

In the next section the non-gaussian distribution of velocity will be related to clusteri- 
zation with the help of a mean field model of driven granular gas. 
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V. A MODEL FOR THE CLUSTERING AND THE NON-GAUSSIAN BEHAVIOR 



In order to shed some light on the relationship between the spatial clusterization and 
the anomalous velocity distribution observed above, we present a simple theoretical model. 
For sake of notation simplicity, we discuss only the Id case. Let us treat the collisions 
in a mean-field like fashion and modify the Langevin dynamics plus collision rules by the 
following set of coupled equations for the velocities: 

JV 



IF—T + Tv^-W^ (47) 



where the second term in the r.h.s. determines the velocity change of the particle i due to 
the collisions with the remaining particles and is chosen to mimick the inelastic behavior. 
This requirement poses some constraints about the form of the function g(v — v')\ 

• The momentum conservation dictates the antisymmetric property, g(v — v') = —g(v' — 
v) 

• The inelasticity of the collision process requires g{y — v')(v — v') < 
The Fokker-Planck equation of (f4?D is 

a 



N rx N 



d t P N (v x , ...,v N ,t) - —(viP N (vi, ...,v N ,t) - — ^ -^Pn{vi, v N , t) + 



8=1 ' 1=1 

N 



N ^ 

Ed 

i=i 



-Vj)PN{vi,...,v N ,t) 



N 



(48) 



^From the above equation, using the fact that in the limit N — > oo the mean field ap- 
proximation holds, one can obtain an evolution equation for the 1-body velocity probability 
distribution which reads: 



J dv'P{v, t)P(v', t)g(v- v') = (49) 



dP(v,t) _ l d{vP{v,t)) _ T F d 2 P{v,t) d_ 
dt t dv t dv 2 dv 

i.e. a sort of self consistent Boltzmann equation. ^From eq. ( fl9"D one observes that the 
quantity: 

dv'P(v', t)g(v - v') = G{v) = (50) 

which is a function of v and a functional of P(v), can be considered as an effective force 
acting on the particle generated by an effective potential U. Integrating once with respect 
to the velocity the stationary version of eq. (1491) one can obtain the following equation: 



V - + ^ + G(v))P(v) = (51) 
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T ( V 2 



The solution of eq. ( [51]) is: 

P(v) oc exp (-— ^- + !/(«) 1 1 (52) 

In order to make some progress we consider the qualitative shape of g(vi — Vj). In eq. 
(f|7|) the effect of collisions between the particles i and j in the unit of time is given by: 

d 2 

j t {vi - Vj)\ coll = —g{vi - vj) (53) 

The variation of momentum in an interval dt can be rewritten as 

6(vi - Vj)\ co a = 5qij ■ (vi - Vj) (54) 

where bqij is the analogue of q = 1 — r in the model discussed in sections and |T|. The 
important difference is that here Sqij represents the effect of all the collisions during dt, and 
thus can be associated to an effective restitution coefficient. Eq. (|5"4]) may be rewritten as: 



- Vj)\ con = XijQ ■ (vi - Vj) (55) 

where Xij is the number of collisions between the i — th and the j — th particles in the unit 
of time. Upon comparing eqns. (|53|) and (|5f|) one obtains an expression for g(vi — Vj): 

g(vi ~ Vj) = - Vj) (56) 

Now it is easy to understand that Xij is a decreasing function of \v{ — Vj\\ indeed, a great 
number of collisions occurs when the pair i, j belongs to a cluster (where \vi — Vj\ is small), 
whereas the two particles rarely collide when they are out of a cluster (and — Vj\ is high). 
We can, therefore, make a rough estimate of g{vi — Vj), that is: 

\v ' — V ' 

\g(vi — Vj)\ ~ — — inside clusters 

\g(vi — Vj) \ ~ \vi — Vjl 13 outside clusters (57) 

where (3' < 1. ^From eq. fl5Tf) it appears that G(v) ~ g(v — v')\ v i =0 as the integration has 



to be performed with respect to the measure P(v',t)dv' that is strongly peaked at v' = 0. 
Finally, one can conclude from the same eq. (B^) that 



v 2 



U(v) v < ^/T g (58) 



r, 



c 



U(v) ~ v p v>JT t 



where (3 — (3' + 1 < 2. It is clear now, looking at eq. (52), that when r < r c (i.e., in the 
non-clusterized regime) the argument of the exponential is dominated by v 2 /r and therefore 
a gaussian is expected for P(v) with variance Tp. In the opposite regime, when r > r c the 
distribution is a gaussian with variance — Tp at low velocities, a simple exponential (if j3 = 1) 
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at high velocities, and a gaussian with variance Tp at extremely high velocities, but this 
very far tails practically cannot be observable. In figure ( ^5|) the tails of the distributions 
of velocities (from the simulation of the model of section [TTI| ) for three different choices of 
parameters are presented: in case (a), when r < r c , we observe a Gaussian distribution; 
in case (b), when r > r c , we can fit the tail with the function exp(— v 3 ' 2 /b), and this is in 
agreement with the analytical calculation performed by van Noije and Ernst ||22|| ; finally in 
case (c), when r >> r c we observe a simple exponential tail, as we may expect from the 
argument presented above. 



VI. CONCLUSIONS AND OPEN PROBLEMS 

In this paper a class of models of granular gases in one and two dimensions has been 
studied by mean of computer simulations and analytical investigations. We think to this 
class of models as the natural, and more physical, extension of previous models in the do- 



main of granular kinetics [16| - [20|. In the models here proposed, by effect of balance 
between Brownian driving and inelastic collisions, one has a good thermodynamic limit; 
furthermore, these models present a rich phenomenology as several regimes are observed by 
tuning the physical parameters, that is the time of viscous interaction r and the coefficient 
of restitution for inelastic collisions r. The two extreme behaviors of those models are the 
Gaussian/homogeneous regime and the non-Gaussian/clusterized one. In the homogeneous 
phase, the system may be described almost as a perfect gas in equilibrium at a temperature 
close to that of the external driving (or a bit lower), showing the absence of densities instabil- 
ities and a Maxwellian distribution of velocities. The out-of-equilibrium phase, on the other 
side, presents strong fluctuations of density (clusters and collapse) with self-similar density 
distribution and a stationary fractal dimension, while there is a strong enhancement of high 
energy tails in the distribution of velocities. This dramatic breaking of the equipartition law 
has to be kept into account in modeling the hydrodynamics of granular media. Furthermore, 
we explained the origin of the different degrees of clusterization by means of a class of balls- 
in-the-boxes models, showing that the effect of inelasticity may be viewed as a bias to the 
transition rates of these random processes: in this context we showed that the non-Gaussian 
distribution of velocities is recovered assuming a local equilibrium with a temperature which 
depends on the local density. The non-Gaussian behavior has been also analytically inves- 
tigated with the help of a model in which the effect of collisions is treated as a mean-field 
force on each particle and using the fact that this force has a different dependence on the 
impact velocity whether the particle is in a cluster or outside of it. Diffusion of particles has 
been also investigated in the simulation of models of section |TJ: no anomalous diffusion has 
been observed. The diffusion coefficients for the non-Gaussian regime have been reported in 
figures (|T6|) and (|17|) showing a weak (and apparently non-monotonic) dependency on the 
restitution coefficient r. A measure of velocity correlation function < v(t)v(t + r) >, which 
appeared not to be a trivial exponential but likely a superposition of different exponentials 
(therefore still integrable in time), has convinced us that, even in the clusterized regime, 
the particles forget their previous velocities rather quickly due to collisions, that is: they 
enter and exit a cluster frequently enough to not affect average diffusion; however, in the 
clusterized regime the diffusion process is dominated by inter-particles collisions, whereas in 
the homogeneous one the diffusion is dominated by the Brownian motion imposed by the 
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model. This is only a rough picture, to be furtherly investigated. 

An important task to accomplish should be the research of an equation of state for this 
class of gases, useful in an eventual hydrodynamic description of them. The observed relation 
between local temperature and density (see figures fllOf) , fl23|) and ( [24]) and the discussion in 
section |TVp is the starting point in this project. Analytical expressions of the pressure have 
to include the usual streaming term p < v 2 > (where p is the local density and < > is an 
ensemble average) as well as a collisional term which is important in the regions where the 
density is high: the streaming term, as a consequence of the scaling < v 2 >~ p" 13 , see the 
first of eqs. P5|), is expected to be proportional to ~ p 1 ' 13 if the picture of local Gaussian 
equilibrium is confirmed ||35|| . 
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FIGURES 
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FIG. 1. Kinetic energy and dissipated energy per particle (as denned in equations (|j) and (|5|)) 
N, in Du et al. [16] model. Particle density N/L = 100 is kept constants and r = 0.99. 
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FIG. 2. The average granular temperature T g and the average dissipated energy per particle 
vs the restitution coefficient r for different values of t and N = 200. From top to bottom: 
= 0.01, r = 2, r = 100, r = 1000. 
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Dissipated energy per particle W versus the granular temperature T g , 



t = 0.01 (a), t = 2 (b), t = 100 (c), r 



for different r 
1000 (d). The dashed lines represent the 
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FIG. 4. Instantaneous density profiles p(x) in two regimes: (a) quasi-equilibrium regime 
(t = 0.01, r = 0.99) and (b) non-equilibrium regime with clusters (r = 100, r = 0.6). In both 
histograms N = 200 and the dashed horizontal lines represent the average density, equal to four 
particles per bin 
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FIG. 5. C(R) vs. R for r = 100, r = 0.6 (top) and r = 100, r = 0.99 (bottom) with N = 200. 
Correlation dimension takes respectively the values d,2 = 0.59 and cfe = 1- 
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FIG. 6. The correlation dimension d,2 vs. r for different values of r: from top to bottom r is 
0.01, 2, 100 and 1000. N = 200 in all simulations 
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FIG. 7. H M /H* M vs. r for different r: form top to bottom r is 0.01, 2, 100 and 1000, with 
N=200 and M=80 {H* M 63) 
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FIG. 8. The rescaled velocity distribution P{v/a) vs. v/a. Pluses (+) are data from simulation 
with t = 100, r = 0.7. Crosses are data with r = 0.01,r = 0.99. The dot-dashed line represents 




the gaussian distribution, while the dashed line represents the fit discussed in section IV 
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FIG. 9. Density distribution vs. m for two different choices of parameters: r = 0.01, r = 0.99 
and r = 100, r = 0.7, N = 300 and M = 100. The former is fitted by a Poisson distribution with 
A = 3, while the latter is fitted by ± e -°- u * m 
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FIG. 10. Box granular temperature TM{m) against m, when r = 0.01, r = 0.99 and r = 100, 
r = 0.7 (in this case th fit m~ ' 5 is plotted). 
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FIG. 11. Snapshot of particle distribution in 2 dimensions in the clusterized regime. ./V = 5000, 
Tc = 0.5, I = 0.63, r = 100, r = 0.01. 
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FIG. 12. C(R) vs. R in the one-dimensional system. TV = 2000, r c = 0.5, / = 0.4. The top 
curve is for r = 100, r = 0.5, the bottom one is for r = 0.01, r = 0.99. Correlation dimension is, 
respectively, 6,2 = 0.55 and d.2 = 1. 
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FIG. 13. C(R) vs. R in 2-d. N = 5000, r c = 0.5, / = 0.71. Top curve is for r = 100, r = 0.5, 
while the bottom one is for r = 0.01, r = 0.99. The correlation dimension is d.2 = 1.45 and c?2 = 2, 
respectively 
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FIG. 14. The granular temperature T g , fractal dimension c?2 and collision rate x vs number of 
particle N, for the model in one dimension, with r c = 0.5, I = 0.4, r = 100 and r = 0.5. 
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FIG. 15. The granular temperature T g , fractal dimension d,2 and collision rate x against number 
of particle N, for the model in two dimensions, with r c = 0.5, / = 0.63, r = 100, r = 0.5 
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FIG. 16. The granular temperature T g , fractal dimension cfe, collision rate \ an d diffusion 
coefficient D against restitution coefficient r, for the model in one dimension, with N = 4000, 
Tc = 0.5, / = 0.4, r = 100 and r = 0.5. 
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FIG. 17. The granular temperature T g , fractal dimension cfe, collision rate x an d diffusion 
coefficient D against restitution coefficient r, for the model in two dimensions, with N = 3000, 
Tc = 0.5, / = 0.63, r = 100, r = 0.5 
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FIG. 18. Distribution of velocities in a gaussian (r = 0.01, r = 0.99) and a non-gaussian regime 
(r = 100, r = 0.5) for the one dimensional system. In both cases N = 2000, r c = 0.5, / = 0.4. The 
dashed curve represents the gaussian. 
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FIG. 19. Distribution of velocities in a gaussian (r = 0.01, r = 0.99) and a non-gaussian regime 
(r = 100, r = 0.5) for the two dimensions case. It is always N = 10000, r c = 0.05, / = 0.22 and 
the dashed line represents the gaussian 
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FIG. 20. Distribution of velocities restricted to number density m = 1 (pluses '+') and m = 5 
(crosses 'X'), in a 2-dimensional case, with N = 10000, r c = 0.05, I = 0.22, r = 100, r = 0.5. 
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FIG. 21. Distribution density /m(^) vs. m for two one-dimensional cases: r = 100, r = 0.5 
and r = 0.01, r = 0.99. In both cases: N = 500, r c = 0.5, Z = 0.4, M = 12000. There are two 
curves superimposed: a Poisson function (with A = N/M ~ 0.04) and m -L95 exp(— 0.26 * m) fit 
for the clusterized regime. 
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FIG. 22. Distribution density /m(^) vs. m for two bidimensional cases: r = 100, r = 0.5 and 
T = 0.01, r = 0.99. In both cases: TV = 10000, r c = 0.05, I = 0.22, M = 3200. There are two 
curves superimposed: a Poisson function (with A = N/M = 3.125) and m~°- 5 exp(— 0.097 * m) fit 
for the clusterized regime. 
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FIG. 23. Box granular temperature TM(m) vs. m for two one-dimensional cases: r = 100, 
r = 0.5 and r = 0.01, r = 0.99. In both cases: N = 500, r c = 0.5, Z = 0.4, M = 12000. The 



gaussian case is constant, while the non-gaussian case is fitted by ~ m 
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FIG. 24. Box granular temperature TM{m) vs. m for two bidimensional cases: r = 100, r = 0.5 
and r = 0.01, r = 0.99. In both cases: N = 10000, r c = 0.05, / = 0.22, M = 3200. The gaussian 
case is constant, while the non-gaussian one is fitted by ~ m 
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FIG. 25. Rescaled distributions of velocities (particular) for three different choices of param- 
eters, in two dimensions: (a) N = 10000, r = 0.01, r = 0.99, with Gaussian fit; (b) N = 3000, 
T = 5, r = 0.5 with the fit ~ exp(-u 3 / 2 /L25); ( c ) N = 10000, r = 100, r = 0.2 with the fit 
~ exp(— v/0.7). In the cases (a) and (c): r c = 0.05, / = 0.22. In the case (b): r c = 0.5 and I = 0.63. 
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